% SizeY,SizeX - Number of Points in X/Y Axis
% phi - Direction of the Gaussian Beam
% res - Resolution in Fractions in Units of lambda
% w0 - Beamwaist in Units of lambda
% xOffset,yOffset - Center of the Beam relative to the center of the Graph
function Gaus = GaussianBeam(SizeX,SizeY,SizeZ,phi,theta,res,w0,xOffset,yOffset)
% Grundvariablen
% Elektrisches Feld
phi = phi*pi/180;
sizeX = 500;
sizeY = 500;
E = zeros(sizeY,sizeX);
R = zeros(sizeY,sizeX);
lambda = 1;
z0 = pi * w0^2/lambda;
k = 2*pi/lambda;
xOffset = xOffset / res;
yOffset = yOffset / res;
for n = 1:sizeY
	y = (n-sizeY/2)/res+0.005;
	for m = 1:sizeX
		x = (m-sizeX/2)/res;
		for o = 1:sizeZ

			r = (cos(phi)*(x-xOffset))-(sin(phi)*(y-yOffset));
			z = (sin(phi)*(x-xOffset))+(cos(phi)*(y-yOffset));
		
			wz = w0 * sqrt(1+(z/z0)^2);
			Rz = z*(1+(z0/z)^2);
			gouy = atan(z/z0);
			E(n,m) = w0/wz * exp(- (r/wz)^2) * exp(-i*k*r^2/(2*Rz))*exp(i*(gouy-k*z));
			%E(n,m) = exp(-i*k*z)*w0/wz*exp(-r;
			Gaus(n,m) = real(E(n,m));
		end	
	end
end
h=surf(Gaus);
set(h, 'edgecolor','none');
